Melting in large sodium clusters: An orbital-free molecular dynamics study. 
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The melting-like transition in sodium clusters Najv, with N=55, 92, and 142 is studied by using 
constant-energy molecular dynamics simulations. An orbital-free version of the Car-Parrinello tech- 
nique is used which scales linearly with system size allowing investigation of the thermal behaviour 
of large clusters. The ground state isomer of Nai42 (an uncomplete three-shell icosahedron) melts 
in two steps: the first one (at « 240 K) is characterized by the high mobility of the atoms located 
on the cluster surface; the second, homogeneous melting (at « 270 K), involves diffusive motion 
of all the atoms across the cluster. For the case of Nag2, the icosahedral structure has a larger 
number of surface vacancies, and melts in two well separated steps, surface melting at ~ 130 K and 
homogeneous melting at « 240 K. Nass, a complete two-shell icosahedron, melts in a single stage 
at ~ 190 K. Our results on homogeneous melting for Nai42 and Na92 are in excellent agreement 
with recent experimental determinations of melting temperatures and latent heats. However, the 
experimentally observed enhancement of the melting temperature around N=55 is not reproduced 
by the calculations. 
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I. INTRODUCTION 

The melting-like transition in finite clusters consisting 
of a small number of atoms is of fundamental interest 
as clusters are often produced in a disordered "liquid" 
state,EJ and it is also relevant to applications of clusters. 
For example, the catalytic activity of small platinum clufe 
ters depends critically on their melting temperatures. El 
Recent experimental advances reveal some details of the 
melting-like transition but, at the. same time, show new 
and interesting features. Martirfl determined the clus- 
ter size dependence of the melting temperature T^ of 
large sodium clusters, composed of thousands of atoms, 
by observing the vanishing of the atomic shell structure 
in the mass spectra upon heating. It was concluded that 
Tm grows with cluster size, but the results did HjOt ex- 
trapolate yet to the T™ of the bulk. Peters et ala per- 
formed X-ray diffraction experiments on large (50 nm) Pb 
clusters and observed the occurrence of surface nielting 
before homogeneous melting. Electron diffrajctionH may 
also help in detecting a surface melting stage.cl Haberland 
and coworkergj have studied the variation with temper- 
ature of the photofragmentation spectra of Najv (N=50- 
200) clusters, and have deduced the melting tempera- 
tures of the clusters. They find that for some cluster 
sizes the melting temperature is a local maximum not 
in exact correspondence with either electronic or atomic 
shell closing numbers, but bracketed by the two, suggest- 
ing that both effects are relevant to the melting process. 

A number of computer simulations of melting in small 
metallic and nonmetallic clusters has been reported, 
the majority of which- employed phenomenological in- 
teratomic potentials .LrQ The use of such parameterized 
potentials allows the consider|a±ion of long dynamical 
trajectories for large clusters.I'Ll Ab initio methods,tHl 
which have also been used, accurately treat the electronic 
structure of the cluster, but are much more expensive 



computationally and are usually restricted to thcj-study 
of small clusters for short dynamical trajectories.Eil Re- 
cently, Rytkonen et a/J13 have performed ah initio molec- 
ular dynamics (aiMD) simulations of the melting of a 
sodium cluster with 40 atoms, but such a large cluster 
required the use of a fast heating rate. [Xhese aiMD 
treatments use the Kohn-Sham (KS) forirO of density 
functional theory (DFT), and orthogonalization of the 
one-electron KS orbitals is the limiting step in their per- 
formance. However, DFT shows that the total energy 
of the electronic system ciasi be expressed in terms of 
just the electronic densityjl^ and orbital-free (OF) ver- 
sions of the aiMD technique based on the electron den- 
sity htjyfi been develappi and employed, both in solid 
statell3ll3 and clusterElTcJ applications. These OF meth- 
ods scale linearly with the system size allowing the study 
of larger clusters for longer simulation times than typ- 
ical aiMD simulations. However, quantum shell effects 
are neglected, so that features associated with electronic 
shell closings jUre not reproduced. 

Previously,E3 we have used the orbital-free molecular 
dynamics (OFMD) method to study the melting pro- 
cess in small sodium clusters. Nag and Na2o, clusters 
outside the range cowered by Haberland's photofragmen- 
tation experiments. eI Here, we report constant energy 
OFMD simulations in a study of the melting-like transi- 
tion in larger clusters, Nass, Na92 and Nai42, which are 
within the size range covered in those experiments, and 
for which a full ah initio treatment of their thermal prop- 
erties would be impractical. Even for the OFMD method 
those large clusters represent a very substantial compu- 
tational effort. The aim of our work is to study the mech- 
anisms by which the melting-like transition proceeds in 
these large clusters. In the next section we briefiy present 
some technical details of the method. The results are pre- 
sented and discussed in section HI and, finally, section IV 
summarizes our main conclusions. 



II. THEORY 

The orbital-free molecular dvjiamics method is a Car- 
Parrinello total energy schemetEl which uses an explicit 
kinetic-energy functional of the electron density, and has 
the electron density as the dynamical variable, as op- 
posed to the KS single particle wavefunctions. In con- 
trast to simulations which use empirical interatomic po- 
tentials, the detailed electronic structure and the elec- 
tronic contribution to the energy and the forces on the 
ions are recalculated efficiently every atomic time-step. 
The main features of the energy functional and the calcu- 
latiorial[S£lieme have been described at length in previous 
work,ll30'Ej'Ej andjietails of our method are as described 
by Aguado et al.EJ In brief, the electronic kinetic en- 
ergy functional of the electron density, n(r), corresponds 
to the gradient exparusion. sucpund the homogeneous limit 
through second ordeitSErEl 



where the first term is the Thomas-Fermi functional 
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and the second is the lowest order gradient correction, 
where T^, the von Weizsacker term, is given by 



T^fd = - 
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The first step of the simulations was the determination 
of low temperature isomers for each of the three cluster 
sizes. For such large clusters it is very difficult to find 
the global minimum because the number of different lo- 
cal minima increases exponentially with the number of 
atoms in the cluster. Instead, one has to adopt struc- 
tures that are likely to have the main characteristics of 
the ground state. We have assumed icosahedral growth. 
Thus, for Nai42 we removed five atoms from the surface of 
a 147 atom three-shell perfect icosahedron. For Nag2, we 
constructed an icosahedral isomer by following the gro" 
ing sequence described by Montejano-Carrizales et al,\ 
and for Nass we took a perfect two-shell icosahedron. We 
have also used dynamical simulated annealing,ll3 to gen- 
erate low temperature isomers, but this procedure always 
led to amorphous structures for Nag2 and Nai42, and to 
a nearly icosahedral structure for Nass . 

Several molecular dynamics simulation runs at differ- 
ent constant energies were performed in order to obtain 
the caloric curve for each icosahedral isomer. The ini- 
tial positions of the atoms for the first run were taken 
by slightly deforming the equilibrium low temperature 
geometry of the isomer. The final configuration of each 
run served as the starting geometry for the next run at a 
different energy. The initial velocities for every new run 
were obtained by scaling the final velocities of the pre- 
ceding run. The total simulation times varied between 8 
and 18 ps for each run at constant energy. 

A number of indicators to locate the melting-like tran- 
sition were employed. Namely, the specific heat defined 
byEl 



The local density approximation is used for exchange and 
correlation.u'Ea In the external field acting on the elec- 
trons, Vexti^) = 'l2n^i^~ ^n)' '^^ " ^O" 

cal pseudopotential of Fiolhais et aL,c3 which reproduces 
well the properties of bulk sodium and has beea shown 
to have good transferability to sodium clusters .E3 

The cluster is placed in a unit cell of a cubic superlat- 
tice, and the set of plane waves periodic in the superlat- 
tice is used as a basis set to emand the valence density. 
Following Car and ParrinelloJliJ the coefficients of that 
expansion are regarded as generalized coordinates of a 
set of fictitious classical particles, and the corresponding 
Lagrange equations of motion for the ions and the elec- 
tron density distribution are solved as described in Ref. 

m 

The calculations for Nag2 and Nai42 used a supercell 
of edge 71 a.u. and the energy cut-off in the plane wave 
expansion of the density was 8 Ryd. For Na55, the cell 
edge was 64 a.u. and the energy cut-off 10 Ryd. In. 
all cases, a 64x64x64 grid was used. Previous testsEHl 
indicate that the cut-offs used give good convergence of 
bond lengths and binding energies. The fictitious mass 
associated with the electron density coefficients ranged 
between 1.0x10^ and 3.3x10* a.u., and the equations of 
motion were integrated using the Verlet algorithmn3 for 
both electrons and ions with a time step ranging from 
At = 0.73 x 10~^^ sec. for the simulations performed 
at the lowest temperatures, to At = 0.34 x 10^^'°^ sec. 
for those at the highest ones. These choices resulted in a 
conservation of the total energy better than 0.1 %. 
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where N is the number of atoms and <>t indicates the 
average along a trajectory; the diffusion coefficient, 

1 d 
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which is obtained from the long time behaviour 
of the mean square displacement < r^it) >~ 

Tk;T,"LiT,f=i[^ii'to, +t) - Ri{ta,)Y, where iit is the 
number of time origins, ^ considered along a trajectory; 
the time evolution of the distance between each atom and 
the center of mass of the cluster 



n{t) =1 R,{t) - Rc,n{t) I 
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and finally, the radial atomic density, averaged over a 
whole dynamical trajectory. 



p{r) = 



dNat{r)/dr 



(7) 



where dNat (r) is the number of atoms at distances from 
the center of mass between r and r + dr. 



III. RESULTS 

The lowest energy structure of sodium clusters of 
medium size is not known. DFT calculations for 
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Na55 performed by Kiimmel et al. using an approx- 
imate structmal model (CAPS model, where the to- 
tal pseudopotepJjial of the ionic skeleton is cylindri- 
cally averaged)t3 give a structure close to icosahedral. 
Near-threshold photoionization mass spectrometry ex- 
periments suggest icosahedral structuresJor large sodium 
clusters with more than 1400 atoms so incomplete 
icosahedral structures are plausible candidates for Na92 
and Nai42. For this reason we have adopted for Nai42 
an isomer obtained by removing five atoms from a per- 
fect three-shell icosahedron. The icosahedral growing se- 
quence in nickeljClusters has been studied by Montejano- 
Carrizales et a?,E3 who have shown that the 12 vertices 
of the outermost shell are the last sites to be occupied. 
Assuming the same growth sequence for sodium clusters, 
we have removed five atoms from the vertex positions of 
Nai47, testing all possibilities, and have then relaxed the 
resulting structures. In the most stable structure thus 
formed, the five vacancies form a pentagon. For Na92, we 
have adopted thc_iimbrella growing model of Montejano- 
Carrizales et aZ.Ej The resulting structure corresponds 
to three complete umbrellas capping a Nass icosahedron. 
Low-temperature dynamical trajectories verify that these 
structures are indeed stable isomers of Na92 and Nai42 . 
The icosahedral isomers are more stable than the lowest 
energy amorphous isomers which were found by simu- 
lated annealing (0.017 eV/atom and 0.020 eV/atom for 
Nag2 and Nai42 respectively). Calvo and Spiegelmann 
have studied sodium clusters in the same size range, utt. 
ing pair potential and tight-binding (TB) calculations, □ 
and have also predicted icosahedral structures for Nass, 
Nags, Nai39 and Nai47- 

For each icosahedral cluster we have calculated the so- 
called caloric curve which is the internal temperature as 
a function of the total energy, where the internal tem- 
perature, is defined as the average of the ionic kinetic 
energyc2l. A thermal phase transition is indicated in the 
caloric curve by a change of slope, the slope being the 
specific heat; the height of the step gives an estimate of 
the latent heat of fusion. However, melting processes are 
more easily recognised as peaks in the specific heat as 
a function of temperature. These have been calculated 
directly from eq. (4) and are shown together with the 
caloric curves in figures 1-3. The specific heat peaks oc- 
cur at the same temperatures as the slope changes of the 
caloric curve giving us confidence in the convergence of 
our results as the two quantities have been obtained in 
different ways. 

The specific heat curve for Nai42 (fig. 1) displays two 
main peaks at Ts « 240 K and T„i « 270 K, suggesting 
a two step melting process, but these are so close that 
only one slope change in the caloric curve can be dis- 
tinguished. Our analysis below shows that homogeneous 
melting occurs at T^ « 270 K in excellent agreement 
with the experiments of Haberland and coworkers,l3 who 
give an approximate value of 280 K. The latent heat of 
fusion estimated from the step at T^ in the caloric curve 
is q,„ ~ 15 meV/atom, again in good agreement with 
the experimental value of ~ 14 meV/atom. However, 
the premelting stage at T=T<; is not detected in the ex- 
periments, but our results are not inconsistent with this 



because the two calculated peaks in the specific heat are 
close together and the height of the first peak is much 
smaller than that of the second; consequently they could 
be difficult to distinguish experimentally. Calvo and 
SpiegelmannEI have performed Monte Carlo (MC) sim- 
ulations using a semiempirical many-atom potential, and 
the lowest-energy isomer they found for Nai39 was also 
an incomplete three-shell icosahedron, in this case with 
8 surface vacancies. They also report two close peaks 
in the specific heat curve indicating a two-step melting 
process, with T^ « 210 K and T™ « 230 K. They con- 
cluded that these two temperatures become closer as the 
cluster size increases, so that for clusters with more than 
about 100 atoms there is effectively just one peak in the 
specific heat and a single-step melting. Tight binding 
(TB) molecular dypUamics calculations were performed by 
the same authorsJj Although the melting temperatures 
were found to be different from those obtained with the 
semiempirical potentials (TB tends to overestimate the 
experimental values, while empirical potentials tend to 
underestimate them), the qualitative picture of melting 
in two close steps was the same. 

The results for Na92 are shown in fig. 2. Two-step 
melting is again observed, with a small prepeak in the 
specific heat at Tg « 130 K and a large peak, correspond- 
ing to homogeneous melting, at Tm ~ 240 K. In this case 
Ts and T,„ are well separated, but the first peak is much 
smaller than the second, which could again account for 
the absence of the prepeak in the experiments. The calcu- 
lated temperature and latent heat for the homogeneous 
melting stage, T™ « 240 K and q™ « 8 meV/atom, 
are again in excellent agreement with the experimental 
values 250 i£ and 7 meV/atom respectively. Calvo and 
Spiegelmanncl arrive at similar conclusions based on MC 
simulations for Nags using phenomenological potentials: 
a small bump near 100 K and a main peak near 180 K. 
Their TB simulations give values for those two tempera- 
tures roughly 100 K higher. 

The experimentsQ indicate a substantial enhancement 
of the melting temperature in the size region around 
N=55 atoms. The reported melting temperature of Na^g 
is 325 K, surprisingly higher than that of Na^42, which is 
a local maximum in the size region of the third icosahe- 
dral shell closing. Our simulations do not reproduce this 
enhancement of T^ for Nass and predict that this clus- 
ter melts in a single stage at T^ « 190 Krifig. 3), a re- 
sult found also by Calvo and Spiegelmann.B The OFMD 
method does not account for electronic quantum-shell ef- 
fects, and full KS calculations may be needed to clarify 
this discrepancy, although it is not clear a priori how elec- 
tronic shell effects could shift the value of T„i by such a 
large amount. Of course, another possibility is that the 
icosahedron is-pot the ground state isomer. However, 
Kiimmel et aO have recently found that the experimen- 
tal photoabsorption spectrum of Na^g is best reproduced 
with a slightly oblate isomer which is close to icosahedral. 
We have also investigated a bcc-like growing sequence 
finding that bcc structures are less stable than icosahe- 
dral ones for all cluster sizes studied. Also, we studied 
the melting behavior of a Nass isomer with bcc structure 
and did not find an enhanced melting temperature for it 
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either. In summary, the discrepancy between experiment 
and theory for Na55 deserves further attention. 

Various quantities have been analyzed in order to in- 
vestigate the nature of the transitions at Tg and T^- The 
short-time averages (sta) of the distances between each 
atom and the center of mass of the cluster, < ri{t) >sta, 
have been calculated, and the cluster evolution during the 
trajectories has been followed visually using computer 
graphics. The < ri{t) >sta curves for Nai42 are pre- 
sented in Figs. 4-6 for three representative temperatures. 
At low temperatures (Fig. 4) the values of < r^it) >sta 
are almost independent of time. The movies show that 
the clusters are solid, the atoms merely vibrating around 
their equilibrium positions. Curve crossings are due to 
oscillatory motion and slight structural relaxations rather 
than diffusive motion. At this low temperature quaside- 
generate groups which are characteristic of the symmetry 
can be distinguished: one line near the centre of mass of 
the cluster identifies the central atom (its position does 
not exactly coincide with the center of mass because of 
the location of the five surface vacancies); 12 lines corre- 
spond to the first icosahedral shell; another 42 complete 
the second shell, within which we can distinguish the 12 
vertex atoms from the rest because their distances to the 
centre of mass are slightly larger; finally, 82 lines describe 
the external shell, where again we can distinguish the 7 
vertex atoms from the rest. 

The radial atomic density distributions with respect 
to the cluster center, p(r), are shown for Nai42 in Fig. 
7. At the lowest temperature, T=30 K, the atoms in the 
icosohedral isomer are distributed in three well separated 
shells, a surface shell and two inner shells; as discussed 
above, subshells form in the second and third shells. The 
shell structure is still present at T=130 K. 

Figures 5 and 7 show that at T=160K the atomic shells 
of the Nai42 cluster are still well defined, but the movies 
reveal isomerization transitions, similar to those found 
at thejpegining of the melting-like transition of Nas and 
Na2o,t2l with no true diffusion. These isomerizations in- 
volve the motion of vacancies in the outer shell, in such a 
way that different isomers are visited which preserve the 
icosahedral structure. The onset of this motion is grad- 
ual and does not lead to features in the specific heat, 
although it is detected in the temperature evolution of 
the diffusion coefficient (see Fig. 8 and discussion be- 
low). The true surface melting stage does not develop in 
the icosahedral isomer until a temperature of Ts w 240 
K is reached. 

Fig. 6 shows the time evolution of < >sta for Nai42 
at a temperature T= 361 K at which the cluster is liq- 
uid with all the atoms diffusing throughout the cluster. 
Some specific cases of atoms that at the begining of the 
simulation are near (far from) the center of mass of the 
cluster and end in a position far from (near) the center 
of mass of the cluster are shown in boldface. The atomic 
density distribution at 280 K, a temperature just above 
the melting point, is nearly uniform across the cluster, a 
radial expansion of the cluster by about 5 bohr units is 
evident, and the surface is more diffuse. 

The < ri{t) >sta curves for Na92 at low temperature 
are qualitatively similar to those of Nai42. Na92 shows 



surface melting at Tg ps 130 K. This temperature is in 
the range where the isomerization processes in Nai42 set 
in, but the larger number of vacancies in the surface shell 
of Na92 allows for more rapid surface diffusion and these 
processes give rise to a distinct peak in the specific heat. 

Nass is a perfect two-shell icosahedron, so surface 
atoms have no empty sites available to move to, and 
diffusion within an atomic shell is almost as difficult as 
diffusion across different shells. When the surface atoms 
have enough energy to exchange positions with one an- 
other they can as easily migrate throughout the whole 
cluster, and melting proceeds in a single stage at 190 K. 
Calvo and SpiegelmanntI have suggested that this one- 
step melting is associated with a large energy gap be- 
tween the ground state icosahedral structure and the 
closest low-lying isomers, but this cannot be a general 
result for perfect icosahedral metallic clusters, as the de- 
tails of melting are material dependent. For example, a 
surface melting sta|ge has been observed in simulations of 
icosahedral Niss.Ea 

The variation of the diffusion coefficient with temper- 
ature is shown in Fig. 8 for Nai42. At temperatures 
less than about 140K, D is close to zero, indicating only 
an oscillatory motion of the atoms. For temperatures 
between 140K and Tg, the diffusion coefficient increases 
indicating that the atoms in the cluster are not undergo- 
ing simple vibrational motion; the atomic motions are, 
nevertheless, of the special kind discussed above that 
preserve the icosahedral structure. The slope of D(T) 
increases appreciably at T^ when surface melting occurs, 
but there is no noticeable feature when the cluster finally 
melts at T^- The features of D(T) for Na92 which are 
not shown here, are very similar: D(T) is very sensitive 
to the surface melting stage, where appreciable diffusive 
motion begins, and the homogeneous melting transition 
is masked by that effect. We conclude that the D(T) 
curve is a good indicator of homogeneous melting only in 
those cases where the surface melting stage is absent, as 
for example in Na55. 

Our results suggest that the melting transition in large 
icosahedral sodium clusters occurs in a smaller tempera- 
ture range than for small clusters such as Nag or Na2o,c2l 
at least near an icosahedral shell closing. Furthermore, 
the size of any prepeak diminishes with respect to the 
main homogenous melting peak as the cluster size in- 
creases, that is as the fraction of atoms that can take 
part in premelting decreases. Consequently, a homoge- 
neous melting temperature can be defined with less am- 
biguity for large clusters. These comments apply to the 
caloric and the specific heat curves, which are the quan- 
tities amenable to experimental measurement. In con- 
trast, microscopic quantities such as the diffusion coeffi- 
cient D or the < ri{t) >sta curves are very sensitive to 
any small reorganization in the atomic arrangement, and 
it is difficult to determine the melting transition from 
the variation of these quantities with temperature. A 
helpful structural, as opposed to thermal, indicator of 
the melting transition in medium-sized or large clusters 
is the shape of the radial atomic density distribution. 
The atomic density displays pronounced shell structure 
at low temperatures which is smoothed out at interme- 
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diate temperatures where the vacancy diffusion and/or 
surface melting mechanisms are present. Above Tm the 
density is flat. 

In figure 9 we compare the calculated values of the 
melting temperature with the exfierimental values. Our 
earlier results for Nag and Na2oc2l are also included, al- 
though for such small sizes there is some ambiguity in 
defining a melting temperature. There is excellent agrees 
ment with the experimental results for Na92 and Nai42.u 
Measurements of the temperature dependence of the pho- 
toabsorption cross sectjeas for Na+ (n=4-16) have re- 
cently been reportedBa-ESl Although the spectra do not 
show evidence of a sharp melting transition, some encour- 
aging comparison between theory and experiment can be 
made. The experimental spectra do not change appre- 
ciably upon increasing the cluster temperature, until at 
T=105 K (the value given as the experimental melting 
temperature of Nas in Fig. 9) the spectra begin to evolve 
in a continuous., way. In our study of the melting be- 
haviour of Nagc^ we found a broad transition starting at 
T=110 K and continuing until T=220 K, at which point 
the "liquid" state was fully developed. This may explain 
the absence of abrupt changes in the experimental pho- 
toabsorption spectrum with temperature. In any case, 
we feel that the good agreement between theory and ex- 
periment may extend to the small sizes. However, our 
method is not expected to give accurate results if oscil- 
lations in the melting temperature with cluster size arise 
as a consequence of electronic shell effects, which is not 
yet known. The discrepancy for Nass remains intriguing. 
In this regard it is noteworthy that our calculated melt- 
ing temperatures for the three large clusters fit precisely 
the expected large N behaviour, Tm(NaAr)=Tm(bulk) -|- 
C/Ns ^ where C is a constant, and yield as a bulk melting 
temperature Tm(bulk)=350 K, which is close to the ob- 
served value of 371 K. A similar extrapolation to the bulk 
melting temperature is not evident in the experimental 
data. 

As it is well known that the specific isomer used to 
start the heating (bpamics can affect the details of the 
melting transition,EZI we have also studied the melting of 
amorphous Na92 and Nai42 clusters that we obtained by 
a simulated annealing. For the amorphous Nai42 cluster 
bulk melting occurred at the same temperature, T=270 
K, as for the icosahedral cluster, while surface melting 
took place at a much lower temperature, T=130 K. How- 
ever, melting of the amorphous Na92 cluster took place 
over a wide temperature, and no sharp transitions were 
detected. We attach little significance to these results as 
the initial structures and the melting behaviour must de- 
pend on the details of the annealing, and the subsequent 
heating. 

IV. DISCUSSION AND CONCLUSIONS 

A few comments regarding the quality of the simula- 
tions are perhaps in order here. The orbital-free represen- 
tation of the atomic interactions, although much more ef- 
ficient than the more accurate KS treatments, is still sub- 
stantially more expensive computationally than a simula- 



tion using phenomenological many-body potentials. Such 
potentials contain a number of parameters that are usu- 
ally chosen by fitting some bulk and/or molecular prop- 
erties. In contrast our model is free of external param- 
eters, although there are approximations in the kinetic 
and exchange-correlation functionals. The orbital-free 
scheme accounts, albeit approximately, for the effects of 
the detailed electronic distribution on the total energy 
and the forces on the ions. We feel that this is partic- 
ularly important in metallic clusters for which a large 
proportion of atoms are on the surface and experience a 
very different electronic environment than an atom in the 
interior. Furthermore, the adjustment of the electronic 
structure and consequently the energy and forces to re- 
arrangements of the ions is also taken into account. But 
the price to be paid for the more accurate description of 
the interactions is a less complete statistical sampling of 
the phase space. The simulation times are substantially 
shorter than those that can be achieved in phenomeno- 
logical simulations. Longer simulations would be needed 
in order to fully converge the heights of the specific heat 
peaks, or in order to observe a van der Waals loop in the 
caloric curves, to mention two examples. But we expect 
that the locations of the various transitions are reliable. 
All the indicators we have used, both thermal and struc- 
tural ones, are in essential agreement on the temperature 
at which the. transitions start. As we discussed in a previ- 
ous paper J£3 longer trajectories may induce just a slight 
lowering in the transition temperatures. 

The melting-like transitions of Nai42, Na92, and 
Nass have been investigated by applying an orbital-free, 
density-functional molecular dynamics method. The 
computational effort which is required is modest in com- 
parison with the traditional Car-Parrinello Molecular 
Dynamics technique based on Kohn-Sham orbitals, that 
would be very costly for clusters of this size. Specifically, 
the computational effort to update the electronic system 
scales linearly with the system size N, in contrast to the 
N'^ scaling of orbital-based methods. This saving allows 
the study of large clusters. However, the price to pay is 
an approximate electron kinetic energy. 

An icosahedral isomer of Nai42 melts in two steps as 
evidenced by the thermal indicators. Nevertheless, there 
are isomerization transitions involving surface defects at 
a temperature as low as 130 K, that preserve the icosa- 
hedral structure and do not give rise to any pronounced 
feature in the caloric curve. The transition at T^ Ki 240 
K from that isomerization regime to a phase in which 
the surface atoms aquire a substantial diffusive motion 
is best described as surface melting. This is followed at 
T,„ « 270 K by homogeneous melting. For Nag2, there 
is a minor peak in at Ts «130K which we associate 
with surface melting. The smaller value of Tg, for this 
cluster compared with Nai42, is due to the less ordered 
surface. Nass, being a perfect two-shell icosahedron with 
no surface defects melts in a single stage at 190 K. In all 
cases, for T>Tm the atoms are able to diffuse through- 
out the cluster volume. Both the calculated T^ at which 
homogeneous melting occurs and the estimated latent 
heat of fusion q^ are in excellent agreement with the ex- 
perimental results of Haberland and coworkers for Nai42 



5 



and Na92 ; our earlier results on the melting of NagcJ are 
also consistent with the variation of the measured opti- 
cal spectrum with temperature. A serious discrepancy 
between theory and experiment remains for Nass. 

We have found that structural quantities obtained 
from the simulations whichj-are very useful in the study 
of melting in small clusters,E2l such as the diffusion coef- 
ficient, are not, in the case of the larger clusters studied 
here, efficient indicators of homogeneous melting, which 
is better located with thermal indicators. A better struc- 
tural indicator is the evolution with temperature of the 
average radial ion density. This quantity flattens when 
homogeneous melting occurs. 



ACKNOWLEDGEMENTS: We would like to 
thank H. Haberland, S. Kiimmel and F. Calvo for 
sending us preprints of their respective works prior 
to publication. This work has been supported by 
DGES (Grants PB95-0720-C02-01 and PB98-D345), 
NATO(Grant CRG.961128) and Junta de Castilla y Leon 
(VA28/99). A. Aguado acknowledges a graduate fellow- 
ship from Junta de Castilla y Leon. M. J. Stott ac- 
knowledges the support of the NSERC of Canada and 
an Iberdrola visiting professorship at the University of 
Valladolid. 

Captions of Figures. 

Figure 1 Caloric and specific heat curves of Nai42, 
taking the internal cluster temperature as the indepen- 
dent variable. The deviation around the mean tempera- 
ture is smaller than the size of the circles. 

Figure 2 Caloric and specific heat curves of Na92 , tak- 
ing the internal cluster temperature as the independent 
variable. The deviation around the mean temperature is 
smaller than the size of the circles. 

Figure 3 Caloric and specific heat curves of Nass , tak- 
ing the internal cluster temperature as the independent 
variable. The deviation around the mean temperature is 
smaller than the size of the circles. 

Figure 4 Short-time averaged distances < ri{t) >sta 
between each atom and the center of mass in Nai42 , as 
functions of time for the icosahedral isomer at T= 30 K. 

Figure 5 Short-time averaged distances < ri{t) >sta 
between each atom and the center of mass in Nai42 , as 
functions of time for the icosahedral isomer at T= 160 K. 
The bold lines follow the evolution of a particular atom 
in the surface shell and another in the outermost core 
shell. 

Figure 6 Short-time averaged distances < ri{t) >sta 
between each atom and the center of mass in Nai42 , as 
functions of time at T= 361 K. The bold lines are to 
guide the eye in following the diffusive behavior of specific 
atoms. 

Figure 7 Time averaged radial atomic densities of the 
icosahedral isomer of Nai42 , at some representative tem- 
peratures. 

Figure 8 Diffusion coefficient as a function of temper- 
ature for the icosahedral isomer of Nai42 . 

Figure 9 Calculated melting temperatures, compared 
with the experimental values. The experimental values 



for the larger cluster sizes are taken from ref. |6|, while 
that for the smallest Nag cluster is taken from ref. ^ 
(see text for details). 
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